using System;
using System.Collections.Generic;
using System.Text;

namespace NumericalRecipes.ModelingOfData
{
  public  class SvdfitMulti
    {

        public SvdfitMulti()
        {
        }
        private double TOL = 1.0e-5;
        private double chisq;
        public double Chi2
        {
            get { return chisq; }
        }
        public void svdfitMulti(double[][] x, double[] y, double[] sig, int ndata, double[] a, int ma,
            double[,] u, double[,] v, double[] w, NumericalRecipes.Delegates.FunctionDoubleArrayToDoubleA funcs)
        {
            NumericalRecipes.SolutionOfLinearAlgebraicEquations.Svbksb obj1
                = new NumericalRecipes.SolutionOfLinearAlgebraicEquations.Svbksb();
            NumericalRecipes.SolutionOfLinearAlgebraicEquations.Svdcmp obj2
                = new NumericalRecipes.SolutionOfLinearAlgebraicEquations.Svdcmp();
            int j, i;
            double wmax, tmp, thresh, sum;
            double[] b = new Double[ndata];
            double[] afunc = new Double[ma];

            for (i = 0; i < ndata; i++)
            {
                afunc = Myfuncs(x[i]);
                tmp = 1.0 / sig[i];
                for (j = 0; j < ma; j++) u[i, j] = afunc[j] * tmp;
                b[i] = y[i] * tmp;
            }
            obj2.svdcmp(u, ndata, ma, w, v);
            wmax = 0.0;
            for (j = 0; j < ma; j++)
                if (w[j] > wmax) wmax = w[j];
            thresh = TOL * wmax;
            for (j = 0; j < ma; j++)
                if (w[j] < thresh) w[j] = 0.0;
            obj1.svbksb(u, w, v, ndata, ma, b, a);
            chisq = 0.0;
            for (i = 0; i < ndata; i++)
            {
                afunc = Myfuncs(x[i]);
                for (sum = 0.0, j = 0; j < ma; j++) sum += a[j] * afunc[j];
                tmp = (y[i] - sum) / sig[i];
                chisq += tmp * tmp;
            }
        }

        // substitutes the x[0..dim-1] values into the function parametrized by
        // a[0..ma-1]
        public double substitute(double[] x, int ndata, double[] a, int ma,
                     NumericalRecipes.Delegates.FunctionDoubleToDoubleA funcs)
        {
            double[] afunc = new Double[ma];
            double sum;
            afunc = Myfuncs(x);
            sum = 0.0;
            for (int j = 0; j < ma; j++) sum += a[j] * afunc[j];
            return sum;
        }
      public double Mysubstitute(double[] x, int ndata, double[] a, int ma)
      {
          double[] afunc = new Double[ma];
          double sum;
          afunc = Myfuncs(x);
          sum = 0.0;
          for (int j = 0; j < ma; j++) sum += a[j] * afunc[j];
          return sum;
      }
        // substitutes the x[0..nData-1, 0..dim-1] values into the function parametrized by
        // a[0..ma-1]
        public double[] substituteVector(double[][] x, int ndata, double[] a, int ma,
                     NumericalRecipes.Delegates.FunctionDoubleArrayToDoubleA funcs)
        {
            double[] afunc = new Double[ma];
            double[] sum = new Double[ndata];
            for (int i = 0; i < ndata; i++)
            {
                afunc = Myfuncs(x[i]);
                sum[i] = 0.0;
                for (int j = 0; j < ma; j++) sum[i] += a[j] * afunc[j];
            }
            return sum;
        }

      public static double[] Myfuncs(double[] x)
      {
          int dimension = x.GetLength(0);
          double[] afunc = new Double[dimension+1];
          afunc[0] = 1;
          for (int i = 1; i <= dimension; i++)
          {
              afunc[i] = x[i-1];

          }        
          return afunc;
      } 
    
    } 
      
}
